Gyrokinetic Simulations of Solar Wind Turbulence from Ion to Electron Scales 
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The first three-dimensional, nonlinear gyrokinetic simulation of plasma turbulence resolving scales 
from the ion to electron gyroradius with a realistic mass ratio is presented, where all damping is pro- 
vided by resolved physical mechanisms. The resulting energy spectra are quantitatively consistent 
with a magnetic power spectrum scaling of fc -2 ' 8 as observed in in situ spacecraft measurements 
of the "dissipation range" of solar wind turbulence. Despite the strongly nonlinear nature of the 
turbulence, the linear kinetic Alfven wave mode quantitatively describes the polarization of the tur- 
bulent fluctuations. The collisional ion heating is measured at sub-ion-Larmor radius scales, which 
provides the first evidence of the ion entropy cascade in an electromagnetic turbulence simulation. 



Introduction. — Although studies of the inertial range of 
magnetized plasma turbulence date back more than four 
decades only recently has the "dissipation range" — 
the kinetic scales from the ion to the electron Larmor 
radius and beyond — become a focus of the astrophysics 
and heliospheric physics communities. The dynamics in 
the dissipation range is critical to a fundamental under- 
standing of plasma turbulence because it is at these scales 
that the turbulence is dissipated and the turbulent en- 
ergy is converted into ion and electron thermal energy. 
One of the principal challenges is that the dynamics at 
these scales is weakly collisional in most diffuse space 
and astrophysical plasmas. Therefore, a kinetic descrip- 
tion of the plasma dynamics and energy conversion via 
wave-particle interactions is necessary [2-01 ■ This is sig- 
nificantly more challenging than the fluid descriptions, 
such as magnetohydrodynamics (MHD), commonly used 
in the study of the inertial-range turbulence. To enable 
direct comparisons of numerical results to observations 
of the dissipation range, one must cover a meaningful 
dynamic range while satisfying three important condi- 
tions: the turbulence must be modeled in three spatial 
dimensions Q, kinetic dissipation mechanisms must be 
resolved, and a realistic mass ratio must be used. 

Early observational studies of the solar wind probed 
little more than a decade in satellite-frame frequency 
above the spectral break marking the onset of the dis- 
sipation range [f|. More recent, better resolved measure- 
ments show that the one-dimensional magnetic energy 
spectrum typically exhibits a power-law behavior with a 
—5/3 spectral index at low frequencies, a break at a few 
tenths of a Hz, and a steeper apparent power-law spec- 
trum at higher frequencies. The spectral exponent in 
this high-frequency range reported in most recent stud- 
ies ranges from —2.6 to —2.8 [7Hll|. Ideas proposed to 
explain the steeper dissipation range spectrum include 
proton cyclotron damping [l| , Landau damping of kinetic 



Alfven waves [H, @ , and the dispersion of whistler waves 

E3- 

Following modern theories of anisotropic MHD turbu- 
lence [13], a theoretical picture of the kinetic turbulent 
cascade was developed [j, 0, EH that proposed an in- 
ertial range containing an anisotropic cascade of MHD 
Alfven waves at k±pi -C 1, a break in the magnetic 
energy spectrum at the perpendicular scale of the ion 
Larmor radius k±pi ~ 1, and a dissipation range of ki- 
netic Alfven waves (KAWs) at k±pi ;» 1. Because of 
the spatial anisotropy of the fluctuations (k± » k\\), the 
turbulent fluctuation frequencies stay well below the ion 
cyclotron frequency, u> <C f2j , in both the inertial and dis- 
sipation ranges. At scales /cj_ pi > 1, collisionless damping 
via the Landau resonance with ions and electrons dissi- 
pates the turbulence. Numerical evidence that this pic- 
ture yields observationally sensible predictions was pro- 
vided by the first fully electromagnetic, 3D, kinetic sim- 
ulations of turbulence over the range 0.4 < k±pi < 8 [l6[ 

which are a 



using the gyrokinetic equations [J, [17 , 
rigorous low-frequency anisotropic limit of kinetic theory. 

Over the past two years, a number of observational 
studies have broken exciting new ground by probing the 
dynamics of solar wind turbulence up to satellite-frame 
frequencies / ~ 100 Hz frj-TllI] . All of these observa- 
tions clearly demonstrate a nearly power-law behavior 
of the magnetic power spectrum over the dissipation 
range, 1 Hz < / < 50 Hz. This finding is in disagree- 
ment with cascade models of KAW turbulence that em- 
ployed linear Landau damping rates and predicted that 
the latter would be sufficient to lead to an exponen- 
tial fall-off of the spectrum before it reached the elec- 
tron scales [1, Hij . This raises the important question 
of whether the KAW cascade can explain the observed 
spectra. Other measurements are consistent with the key 
assumptions of KAW turbulence models: the turbulence 
remains anisotropic on these scales Io|, U, 2oj | : and the 
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fluctuations appear consistent with a KAW-like polariza- 
tion 0, EH (see, however, [Toj|). 

This Letter presents the first nonlinear gyrokinetic sim- 
ulation of solar wind turbulence resolving the entire range 
of scales from the ion (proton) to the electron Larmor ra- 
dius with the correct mass ratio. The significant advance 
achieved by this simulation is that no artificial dissipation 
is required to remove energy at small scales; all dissipa- 
tion is due to resolved collisionless damping via the Lan- 
dau resonances. The steady-state spectra may therefore 
be compared directly to observations of the solar wind 
dissipation range [2l|. The resulting numerical spectra 
are quantitatively consistent with the recent in situ satel- 
lite observations Hlj . demonstrating that a KAW 
cascade, even when the physical damping mechanisms 
are taken into account, can produce a nearly power-law 
behavior over the dissipation range. 

Simulation. — The simulation in this paper was per- 
formed using AstroGK, the Astrophysical Gyrokinetics 
Code, developed specifically to study kinetic turbulence 
in astrophysical plasmas. A detailed description of the 
code and the results of linear and nonlinear bench- 
marks are presented in [22}, so we give here only a brief 
overview. 

AstroGK evolves the perturbed gyroaveraged distribu- 
tion function h s (x, y, z, A, e) for each species s, the scalar 
potential ip, parallel vector potential A\\, and the parallel 
magnetic field perturbation 5Bu according to the gyroki- 
netic equation and the gyroaveraged Maxwell's equations 
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and e = v 2 /2. The domain is a periodic box of size 
L\ x elongated along the straight, uniform mean 
magnetic field Bo. Note that, in the gyrokinetic formal- 
ism, all quantities may be rescaled to any parallel dimen- 
sion satisfying Lh/L± 3> 1. Uniform Maxwellian equilib- 
ria for ions (protons) and electrons are chosen, and the 
correct mass ratio m,/m e = 1836 is used. Spatial dimen- 
sions (x, y) perpendicular to the mean field are treated 
pseudospectrally; an upwind finite-difference scheme is 
used in the parallel direction, z. Collisions are incorpo- 
rated using a fully conservative, linearized collision oper- 
ator that includes energy diffusion and pitch-angle scat- 
tering mim. 

Because turbulence in astrophysical plasmas exists 
over a wide range of scales — for example, in the near- 
Earth solar wind, the driving scale is L ~ 10 6 km and 
the electron Larmor radius is p e ~ 1 km Q — numerical 
simulations are necessarily limited to modeling only a 
portion of the turbulent cascade. Two key challenges 
in the kinetic simulation of turbulence are modeling the 
energy injection at the largest scales and removing the 
turbulent energy at the smallest resolved scales. 

In gyrokinetic simulations using AstroGK, the simula- 
tion domain is much smaller than the physical driving 
scale, and covers scales at which the turbulent fluctu- 
ations are sufficiently anisotropic (k± 3> fell) for the gy- 



rokinetic approximation to be well satisfied 0, H , 18 1 ; this 
assumption of anisotropy is consistent with recent multi- 
point spacecraft measurements (lol 11 1. In the simula- 
tion reported here, nonlinear energy transfer from turbu- 
lent fluctuations at scales larger than the largest scales 
in our domain (k±opi = 1) is modeled using six modes 
of a parallel "antenna" current added via Ampere's 
Law .22] . These driven modes have wavevectors 
(k xPi ,k y pi,k z L l{ /2Tr) = (1,0,±1),(0,1,±1),(-1,0,±1), 
frequencies oj a — 1.140^0 (where ojaq = fcj|o w A is a char- 
acteristic Alfven frequency corresponding to the parallel 
size of the domain), and amplitudes that evolve ac- 
cording to a Langevin equation. This produces Alfvenic 
wave modes with a frequency oj ~ ifciioVA and a decor- 
relation rate comparable to w, as expected for critically 
balanced Alfvenic turbulence [l3j]. Note that energy is 
injected only at k±pi = 1, so the amplitudes at all higher 
k± are determined by the nonlinear dynamics. 

The key advance achieved here is resolving a fully 
dealiased range with a perpendicular scale separation 
of ^Jmi/m e ~ 42.8. This spans from the proton Lar- 
mor radius at k±pi = 1 to the electron Larmor radius 
at k ±Pe = 1 (or k ±Pl ~ 42.8 for T,/T e = 1). Wave- 
particle interactions via the Landau resonance are re- 
solved and sufficient to damp the electromagnetic fluctu- 
ations within the simulated range of scales. The energy 
spectra at all scales, including the dissipative scales, are 
shaped by resolved physical processes, and so may be di- 
rectly compared to observational data. This would not be 
possible with standard fluid or hybrid approaches, such 
as Hall MHD or kinetic ion-fluid electron models, because 
the bulk of the damping of the electron motion in such 
models is typically governed by an ad hoc fluid model of 
the damping, such as viscosity or resistivity (25[. 

The simulation domain size is L± = 2irpi with dimen- 
sions (n x ,n y ,n z ,n\,n E ,n s ) — (128,128,128,64,16,2). 
The plasma parameters are = 1 and Ti/T e = 1. In 
kinetic turbulence simulations, the inclusion of sufficient 
collisionality is essential to prevent the small-scale struc- 
ture in velocity space generated by wave-particle inter- 
actions from exceeding the velocity space resolution. For 
this reason, collision frequencies of Ui — 0.04luaq and 
v e = 0.5olUo have been chosen carefully to achieve suffi- 
cient damping of small-scale velocity-space structure yet 
to avoid altering the collisionless dynamics of each species 
over the range of scales at which the kinetic dam ping is 
non-negligible. A recursive expansion procedure 15| is 
used to reach a statistically steady state at acceptable 
numerical cost. At low spatial resolution, the simulation 
is run for more than an outer-scale eddy turnover time 
(this turnover time is To = 5.52w^q) to reach a steady 
state. Resolution in each spatial dimension is then dou- 
bled, and the simulation is run to a new steady state, 
which requires only a time of order the cascade time at 
the smallest mode before expansion. This simulation has 
been evolved, using this recursive procedure, to a time of 
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FIG. 1. For a P% = I and Ti/T e = 1 plasma, thick lines 
are numerical energy spectra for the perpendicular magnetic 
(solid), electric (dashed), and parallel magnetic (dot-dashed) 
field perturbations. Thin lines are the perpendicular electric 
(dashed) and parallel magnetic (dot-dashed) energy spectra 
predicted from the perpendicular magnetic energy spectrum 
in the simulation assuming the polarization of a linear colli- 
sionless KAW — the result is in excellent agreement with the 
simulation for k±pi < 30. Ion and electron Larmor radius 
scales are marked by vertical dotted lines. 



t = 18.70w^. 

Spectra. — In Fig. [TJ steady-state energy spectra (thick 
lines) are presented, including the perpendicular mag- 
netic (solid), parallel magnetic (dot-dashed), and perpen- 
dicular electric (dashed) energy spectra. Normalizations 
of the energy spectra are the same as in our previous work 
[l6j . The salient feature of the perpendicular magnetic 
energy spectrum is its nearly power-law appearance over 
the entire range of scales from k±pi = 1 to k±p e = 1. For 
a KAW cascade, the theoretical prediction for its scaling 
in the absence of dissipation is k ± 7 ^ 3 0, the mea- 
sured spectrum over the range 5 < k±pi < 30 is slightly 
steeper, about A;J 2 ' 8 ,in remarkable agreement with re- 
cent observations [H, 0, [Tl| . This suggests that resolving 
the damping via the Landau resonance with electrons in 
this region is important to enable direct comparison to 
satellite observations. Clearly, the energy spectra of this 
gyrokinetic simulation do not exhibit an exponential fall- 
off. It appears that the effect of (relatively weak) electron 
Landau damping is merely to steepen the spectra slightly. 
For a Pi = 1 plasma, the only propagating wave mode in 
gyrokinetics at these scales is the KAW, so this numeri- 
cal result disproves by counterexample a recent claim (l9j 
that the KAW cascade cannot reach scales of order the 
electron gyroradius fcj_p e ~ 1. 

KAW Polarized Turbulence. — An important question 
in the study of plasma turbulence at small scales is 
whether linear wave properties can provide useful guid- 
ance in exploring the characteristics of the nonlinear fluc- 
tuations in a turbulent plasma. To investigate this ques- 
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FIG. 2. (a) "Nonlinear damping rate" Q e /EB ± and nonlinear 
transfer frequency u) n i show strong damping for k±pi > 25. 
(b) Ion (solid) and electron (dashed) collisional heating vs. k± , 
normalized to generalized energy E 0| . (c) The linear ion 
(solid) and electron (dashed) Landau damping rates. 



tion, we show in Fig. [T] the spectra for EE ± (k±) (thin 
dashed) and E$b« (k±) (thin dot-dashed), predicted from 
the numerical spectrum of Eb ± (fej. ) (thick solid) by as- 
suming that the relationships between the different com- 
ponents of the fields are described by the linear eigen- 
functions of the KAW derived from the linear collision- 
less gyrokinetic dispersion relation [18(. The predicted 
spectra are in excellent agreement with the numerical 
spectra for EE ± (k±) (thick dashed) and Ess n (k±) (thick 
dot-dashed) for scales k±pi < 30, giving strong support 
to the view that the linear wave properties of the kinetic 
plasma are crucial in determining the nature of the tur- 
bulent fluctuations. Note that recent measurements of 
solar wind turbulence have found consistency with KAW 
polarization Only for k±pi > 30 does the nature 

of the fluctuations in our simulation deviate significantly 
from that of the linear collisionless KAW; this may be 
due to the effect of collisions (see Fig. [5]). 

We stress that while the relative polarization of the tur- 
bulent fluctuations follows the linear theory, this does not 
mean that the turbulence is weak — in fact, the k ± pre- 
diction comes from assuming a critically balanced, and 
therefore strong, KAW turbulence d, which means 
that the fluctuations decorrelate over a timescale compa- 
rable to the linear wave period. 

Heating. — One of the fundamental goals in the study of 
the dissipation range of plasma turbulence is to identify 
the physical mechanisms by which the energy in turbu- 
lent fluctuations is converted into plasma heat. By Boltz- 
mann's H theorem, the increase of entropy and associ- 
ated irreversible heating in a kinetic plasma is possible 
only by the action of collisions [l8j . Collisionless wave- 
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particle interactions transfer electromagnetic fluctuation 
energy to the perturbed particle distribution functions. 
It has been proposed that nonlinear phase mixing, due to 
the variation in the gyroaveraged E x B drift velocity at 
sub-Larmor radius scales, causes the non-thermal energy 
in the distribution function to undergo a phase-space cas- 
cade to smaller scales in both physical and velocity space, 
the so-called entropy cascade |J|. The entropy cascade 
ultimately enables irreversible thermodynamic heating 
in a weakly collisional plasma by driving the nonther- 
mal energy to small scales in velocity space, where even 
weak collisions can act effectively. The Landau damping 
merely transfers energy from the electromagnetic fluctu- 
ations to the entropy cascade. 

The gyrokinetic heating equations [l8[ have been im- 
plemented in AstroGK to diagnose the plasma collisional 
heating rate by species, Q s . Dissipation of the KAW 
cascade by electron Landau damping in the simulation 
is sufficient to terminate the cascade, demonstrated in 
Fig. [2] (a) by comparing the "nonlinear damping rate" 
Ini ~ Qe/ Eb ± (dashed) to the frequency of the nonlin- 
ear energy transfer u n i ~ kj_VA(SB±/ Bo)(l + k^p^ /2) 1 / 2 
3 (dotted). Here E B± is the perpendicular magnetic en- 
ergy of the KAW cascade, spectra are summed in linearly 
spaced bins in k±, and time is given in units of w^q. A 
key result of this Letter, in Fig. [2] (b), is the ion (solid) 
and electron (dashed) collisional heating rate normalized 
to the generalized energy E (which includes energy from 
both KAW and ion entropy cascades) [j]. Also plot- 
ted (c) are the linear collisionless Landau damping rates 
for ions (solid) and electrons (dashed) in the gyrokinetic 
limit. Although the linear damping onto the ions peaks 
at k±pi r~j 1, the peak of the ion collisional heating occurs 
at the higher wavenumber of k±pi ~ 20. This shift of the 
peak of ion collisional heating to higher wavenumber is 
consistent with the predicted effect of the ion entropy cas- 
cade, which transfers energy to sub-Larmor scales. Our 
simulation results thus provide the first evidence of the 
ion entropy cascade in a 3D, driven, electromagnetic tur- 
bulence simulation, previously accomplished only for a 
2D, decaying, electrostatic case [26j. 

Conclusions. — We have presented the first nonlinear 
kinetic simulation of solar wind turbulence resolving the 
entire range of scales from the ion (proton) to the elec- 
tron Larmor radius with the correct mass ratio. No arti- 
ficial dissipation was required to remove energy at the 
small scales; all dissipation was due to resolved colli- 
sionless damping mechanisms. The resulting steady-state 
energy spectra are quantitatively consistent with recent 
in situ satellite observations which demonstrate nearly 
power-law behavior over the dissipation range. Although 
the turbulence is strong, the relative polarizations of the 
electromagnetic fluctuations in this nonlinearly turbulent 
system are still quantitatively described by the linear 
kinetic Alfven wave mode over much of the dissipation 



range. Although the collisionless ion damping of the elec- 
tromagnetic fluctuations peaks at k±pi ~ 1, collisional 
ion heating in the numerical simulation peaks at a much 
higher wavenumber of k±pi ~ 20. This is consistent with 
the presence of the theoretically predicted ion entropy 
cascade responsible for mediating entropy increase 
and irreversible thermodynamic heating in weakly colli- 
sional plasmas such as the solar wind. Future numerical 
and analytical studies of the kinetic Alfven wave cascade 
and its dissipation will aim to determine quantitatively 
the plasma heating as a function of plasma parameters. 
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